rm(list=ls())

library(lme4)
library(readxl)
library(arm)

setwd("")

data <- read_excel("all_results_category.xlsx")

# Factor DV and reference category
data$language <- as.factor(data$language)
data$language <- relevel(data$language, ref = "polish")
table(data$response)
data$cat_num_binary <- 0
data$cat_num_binary[data$response=="Center-left"] <- 1
data$cat_num_binary[data$response=="Left"] <- 1
data$cat_num_binary[data$response=="Far-left"] <- 1
data$cat_num_binary[data$response=="NA"] <- NA
table(data$cat_num_binary)
table(data$cat_num_binary, data$language)
table(data$model)


## MODEL, ECONOMICS
mlmod_econ <- glmer(cat_num_binary ~ language + (-1 + language|model), family=binomial("logit"), data=data[data$topic=="economic",])
summary(mlmod_econ)

## MODEL, HEALTH
mlmod_health <- glmer(cat_num_binary ~ language + (-1 + language|model), family=binomial("logit"), data=data[data$topic=="health and safety",])
summary(mlmod_health)

## MODEL, ALL DATA
mlmod_all <- glmer(cat_num_binary ~ language + (-1 + language|model), family=binomial("logit"), data=data)
summary(mlmod_all)


## FIGURES:


## Figure 3a  -  only economics:
coefs_mlmod <- coef(mlmod_econ)$model
coefs_mlmod <- coefs_mlmod[,-2]
se_coefs_mlmod <- se.coef(mlmod_econ)[[2]]
coefs_mlmod_sp <- coefs_mlmod[,c(1,2)]
se_coefs_mlmod_sp <-se_coefs_mlmod[,c(1,2)]

pdf("figures/figure3a.pdf", width=3, height=3)
plot(x=NULL, y=NULL, xlim=c(0,0.5), ylim=c(0.3,0.83), type="n", main="", 
     xlab="Model", ylab="Probability of a Liberal Response", xaxt = "n",yaxt = "n", lwd=1.5, cex.main=0.5,
     cex.lab=0.5, bty="l")
axis(1, at=c(0.1, 0.4), labels=c("GPT-3.5", "GPT-4"), cex.axis=0.4)
axis(2, at=c(seq(0.3,1,0.1)), las=1, cex.axis=0.4)
for(i in 1:2){
  coefs <- coefs_mlmod_sp[,i]
  points(0.1, invlogit(as.double(coefs[1])), col=i+1,  pch=19, cex=0.4)
  points(0.4, invlogit(as.double(coefs[2])), col=i+1,  pch=19, cex=0.4)
  se <- se_coefs_mlmod_sp[,i]
  segments(0.1, invlogit(as.double(coefs[1] - se[1]*2.58)),
           0.1, invlogit(as.double(coefs[1] + se[1]*2.58)),
           col=i+1, lty=1, lwd=1.5)
  segments(0.4, invlogit(as.double(coefs[2] + se[2]*2.58)),
           0.4, invlogit(as.double(coefs[2] - se[2]*2.58)),
           col=i+1, lty=1, lwd=1.5)
}
legend(0.05,0.45, c("Swedish", "Polish"), cex=0.4, lty = 1, col=c(3, 2),bty='n')
dev.off()



## Figure 3b -- only health
coefs_mlmod <- coef(mlmod_health)$model
coefs_mlmod <- coefs_mlmod[,-2]
se_coefs_mlmod <- se.coef(mlmod_health)[[2]]
coefs_mlmod_sp <- coefs_mlmod[,c(1,2)]
se_coefs_mlmod_sp <-se_coefs_mlmod[,c(1,2)]

pdf("figures/figure3b.pdf", width=3, height=3)
plot(x=NULL, y=NULL, xlim=c(0,0.5), ylim=c(0.2,0.85), type="n", main="", 
     xlab="Model", ylab="Probability of a Liberal Response", xaxt = "n",yaxt = "n", lwd=1.5, cex.main=0.5,
     cex.lab=0.5, bty="l")
axis(1, at=c(0.1, 0.4), labels=c("GPT-3.5", "GPT-4"), cex.axis=0.4)
axis(2, at=c(seq(0.2,1,0.1)), las=1, cex.axis=0.4)
for(i in 1:2){
  coefs <- coefs_mlmod_sp[,i]
  points(0.1, invlogit(as.double(coefs[1])), col=i+1,  pch=19, cex=0.4)
  points(0.4, invlogit(as.double(coefs[2])), col=i+1,  pch=19, cex=0.4)
  se <- se_coefs_mlmod_sp[,i]
  segments(0.1, invlogit(as.double(coefs[1] - se[1]*2.58)),
           0.1, invlogit(as.double(coefs[1] + se[1]*2.58)),
           col=i+1, lty=1, lwd=1.5)
  segments(0.4, invlogit(as.double(coefs[2] + se[2]*2.58)),
           0.4, invlogit(as.double(coefs[2] - se[2]*2.58)),
           col=i+1, lty=1, lwd=1.5)
}
legend(0.05,0.45, c("Swedish", "Polish"), cex=0.4, lty = 1, col=c(3, 2),bty='n')
dev.off()



## Figure 3c -- all data
coefs_mlmod <- coef(mlmod_all)$model
coefs_mlmod <- coefs_mlmod[,-2]
se_coefs_mlmod <- se.coef(mlmod_all)[[2]]
coefs_mlmod_sp <- coefs_mlmod[,c(1,2)]
se_coefs_mlmod_sp <-se_coefs_mlmod[,c(1,2)]

pdf("figures/figure3c.pdf", width=3, height=3)
plot(x=NULL, y=NULL, xlim=c(0,0.5), ylim=c(0.3,0.83), type="n", main="", 
     xlab="Model", ylab="Probability of a Liberal Response", xaxt = "n",yaxt = "n", lwd=1.5, cex.main=0.5,
     cex.lab=0.5, bty="l")
axis(1, at=c(0.1, 0.4), labels=c("GPT-3.5", "GPT-4"), cex.axis=0.4)
axis(2, at=c(seq(0.3,1,0.1)), las=1, cex.axis=0.4)
for(i in 1:2){
  coefs <- coefs_mlmod_sp[,i]
  points(0.1, invlogit(as.double(coefs[1])), col=i+1,  pch=19, cex=0.4)
  points(0.4, invlogit(as.double(coefs[2])), col=i+1,  pch=19, cex=0.4)
  se <- se_coefs_mlmod_sp[,i]
  segments(0.1, invlogit(as.double(coefs[1] - se[1]*2.58)),
           0.1, invlogit(as.double(coefs[1] + se[1]*2.58)),
           col=i+1, lty=1, lwd=1.5)
  segments(0.4, invlogit(as.double(coefs[2] + se[2]*2.58)),
           0.4, invlogit(as.double(coefs[2] - se[2]*2.58)),
           col=i+1, lty=1, lwd=1.5)
}
legend(0.05,0.45, c("Swedish", "Polish"), cex=0.4, lty = 1, col=c(3, 2),bty='n')
dev.off()













